Figure 4
rstudio
September 11, 2024 14:10:32 CEST
Introduction
Load libraries
Define colors
Load single cell RNA-Seq data
Figure 4a
Figure 4b
Figure 4c
Figure 4d
Figure 4e
Redefine scRepertoire function
To be able to change few ploting parameters of the
clonalHomeostasis function, we copy/paste and modified some
function from scRepertoire
Open clonotype data
n= 29690 TCR with EXACTLY one alpha chain CDR3 and one beta chain CDR3
Add transcriptomic data
if(params$accessibility == "unlock"){
srat <- RenameCells(srat,new.names = paste0(srat$patient, "_",gsub("(_1|_2|_3|_4|_5|_6|_7|_8|_9)*", "",colnames(srat))))
srat <- RenameCells(srat,new.names = gsub("_post", "",colnames(srat)))
srat$sample <- gsub("_post", "",srat$patient)
seurat <- combineExpression(combined, srat,
cloneTypes = c(Rare = .00005, Small = .0005,
Medium = .005, Large = .05, Hyperexpanded = 1),
cloneCall = "strict",
group.by = "sample", chain = "both",
proportion = TRUE,
filterNA = TRUE)
s <- subset(seurat, subset = anno_l1 %in% c("CD8+ T cells","CD4+ T cells", "Proliferating cells", "Tregs", "Natural killer cells"))
s$cloneType <- droplevels(s$cloneType)
s@meta.data$Pathological.Response <- factor(s@meta.data$Pathological.Response, levels = c("pCR", "non-pCR"))
} else{
print("Please request access to the BCR-Seq data.")
}We mapped 26528 TCR to the scRNA-Seq data, including 23099 to annotated T cells (10699 CD4 T cells, 1951 Tregs and 9517 CD8 T cells). Only 3 patients do not have hyperexpanded clones (2 non-pCR and 1 pCR).
if(params$accessibility == "unlock"){
ss <- subset(seurat, subset = anno_l1 %in% c("CD8+ T cells") )
ss@meta.data$cloneType <- factor(ss@meta.data$cloneType, levels = c("Hyperexpanded (0.05 < X <= 1)", "Large (0.005 < X <= 0.05)", "Medium (5e-04 < X <= 0.005)", "Small (5e-05 < X <= 5e-04)"))
# define cell group membership
Idents(ss) <- ss$cloneType
de_markers <- DElegate::FindAllMarkers2(ss, replicate_column = "patient", method = "edger", min_fc = log2(1), min_rate = 0.1)
signature_h <- de_markers$feature[de_markers$group1=="Hyperexpanded (0.05 < X <= 1)"]
background <- rownames(ss)
# enricher for cnet plot
GO_H_gene_sets = msigdbr(species = "human", category = "H")
msigdbr_t2g = GO_H_gene_sets %>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()
ego_h <- enricher(gene = signature_h, universe = background, TERM2GENE = msigdbr_t2g)
b <- barplot(ego_h, title = "Hallmarks´enrichment of hyperexpanded clones")
tmp <- b$data
tmp$ID <- str_trunc(tmp$ID, 100, "right")
tmp$ID <- factor(tmp$ID, levels = tmp$ID[order(tmp$Count)])
p <- ggplot(tmp, aes(Count, ID)) +
geom_bar(stat = "identity", color="black", fill = colors_clonotype["Hyperexpanded (0.05 < X <= 1)"]) +
theme_classic()+
geom_text(
aes(label = (paste( "padj=", formatC(p.adjust, format = "e", digits = 2)))),
color = "black",
size = 6,
hjust=1,
position = position_dodge(0.5) ) +
theme(text = element_text(size=24)) +
ggtitle(" Hyperexpanded \n CD8+ T cells clones")
print(p)
DT::datatable(p$data,
caption = ("Figure 4e"),
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c( 'csv', 'excel')))
} else{
print("Please request access to the BCR-Seq data.")
}
Figure 4f
We maped 1286 TRA out of 23099 detected TRA (7.3%) We maped 435 unique TRA out of 10741 unique TRA chains (4%).
if(params$accessibility == "unlock"){
McPAS <- utils::read.csv(params$McPAS)
head(McPAS)
a <- colsplit(s$CTaa, "_", name=c("CTaa_TRA", "CTaa_TRB"))
rownames(a) <- colnames(s)
a$rid <- rownames(a)
c <- McPAS[McPAS$CDR3.alpha.aa %in% na.omit(a$CTaa_TRA), c("CDR3.alpha.aa", "Category", "Pathology")]
head(c)
colnames(c) <- c("CTaa_TRA", "Category", "Pathology")
joint <- left_join(a, c, by = "CTaa_TRA")
joint <- joint[joint$rid %in% colnames(s),]
joint <- joint[!duplicated(joint$rid),]
rownames(joint) <- (joint$rid)
s <- AddMetaData(s, joint)
Idents(s) <- "Pathology"
cells <- WhichCells(s, ident = c("Human herpes virus 1" ) )
Idents(s) <- "Category"
s <- SetIdent(s, cells = cells, value = paste0("hHSV: ", s@meta.data$CTaa[s@meta.data$Pathology %in% c("Human herpes virus 1")]))
s$Category <- Idents(s)
Idents(s) <- "Pathology"
cells <- WhichCells(s, ident = c("Herpes simplex virus 2 (HSV2)") )
Idents(s) <- "Category"
s <- SetIdent(s, cells = cells, value = paste0("hHSV: ", s@meta.data$CTaa[s@meta.data$Pathology %in% c("Herpes simplex virus 2 (HSV2)")]))
s$Category <- Idents(s)
} else{
print("Please request access to the BCR-Seq data.")
}if(params$accessibility == "unlock"){
Idents(s) <- "Category"
cells <- WhichCells(s, ident = c("Pathogens","Autoimmune", "Allergy" ) )
Idents(s) <- "Category"
s <- SetIdent(s, cells = cells, value = "Other")
s$Category <- Idents(s)
s <- subset(s, subset = Category != "Other")
s$Category <- droplevels(s$Category)
p3 <- SCpubr::do_BarPlot(subset(s, subset = Category != "NA" & Category != "Other" & Pathological.Response == "pCR"),
group.by = "Category",
split.by = "patient",
order=FALSE,
legend.position = "right",
colors.use = colors_McPAS,
font.size = 28, flip = TRUE) + xlab("pCR")
p4 <- SCpubr::do_BarPlot(subset(s, subset = Category != "NA" & Category != "Other"& Pathological.Response == "non-pCR"),
group.by = "Category",
split.by = "patient",
order=FALSE,
legend.position = "right",
colors.use = colors_McPAS,
font.size = 28, flip = TRUE) + xlab("non-pCR")
p3$labels$y <- " "
p <- p3 + p4 + plot_layout(ncol = 1, heights = c(1,2), guides = "collect") + ylab("Number of mapped TCR")
print(p)
DT::datatable(rbind(p3$data, p4$data),
caption = ("Figure 4f"),
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c( 'csv', 'excel')))
} else{
print("Please request access to the BCR-Seq data.")
}
Session Info
## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Vienna
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] stringr_1.5.1 msigdbr_7.5.1 DOSE_3.26.2 org.Hs.eg.db_3.17.0
## [5] AnnotationDbi_1.62.2 IRanges_2.34.1 S4Vectors_0.38.2 Biobase_2.60.0
## [9] BiocGenerics_0.46.0 clusterProfiler_4.8.3 enrichplot_1.20.3 scales_1.3.0
## [13] RColorBrewer_1.1-3 ggnewscale_0.4.10 tidyr_1.3.1 scRepertoire_1.10.1
## [17] dittoSeq_1.12.2 canceRbits_0.1.6 ggpubr_0.6.0.999 ggplot2_3.5.1
## [21] viridis_0.6.5 viridisLite_0.4.2 reshape2_1.4.4 tibble_3.2.1
## [25] SCpubr_2.0.2 DT_0.32 patchwork_1.2.0 dplyr_1.1.4
## [29] Seurat_5.0.3 SeuratObject_5.0.1 sp_2.1-3
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.4 matrixStats_1.2.0 spatstat.sparse_3.0-3
## [4] bitops_1.0-7 HDO.db_0.99.1 httr_1.4.7
## [7] doParallel_1.0.17 tools_4.3.0 sctransform_0.4.1
## [10] backports_1.4.1 utf8_1.2.4 R6_2.5.1
## [13] vegan_2.6-4 lazyeval_0.2.2 uwot_0.1.16
## [16] mgcv_1.9-1 permute_0.9-7 withr_3.0.0
## [19] gridExtra_2.3 progressr_0.14.0 cli_3.6.2
## [22] spatstat.explore_3.2-7 fastDummies_1.7.3 scatterpie_0.2.1
## [25] isoband_0.2.7 labeling_0.4.3 sass_0.4.9
## [28] spatstat.data_3.0-4 ggridges_0.5.6 pbapply_1.7-2
## [31] yulab.utils_0.1.4 gson_0.1.0 stringdist_0.9.12
## [34] parallelly_1.37.1 limma_3.56.2 RSQLite_2.3.5
## [37] VGAM_1.1-10 rstudioapi_0.16.0 generics_0.1.3
## [40] gridGraphics_0.5-1 ica_1.0-3 spatstat.random_3.2-3
## [43] crosstalk_1.2.1 car_3.1-2 GO.db_3.17.0
## [46] Matrix_1.6-5 ggbeeswarm_0.7.2 fansi_1.0.6
## [49] abind_1.4-5 lifecycle_1.0.4 edgeR_3.42.4
## [52] yaml_2.3.8 carData_3.0-5 SummarizedExperiment_1.30.2
## [55] qvalue_2.32.0 Rtsne_0.17 blob_1.2.4
## [58] grid_4.3.0 promises_1.2.1 crayon_1.5.2
## [61] miniUI_0.1.1.1 lattice_0.22-6 cowplot_1.1.3
## [64] KEGGREST_1.40.1 pillar_1.9.0 knitr_1.45
## [67] fgsea_1.26.0 GenomicRanges_1.52.1 future.apply_1.11.1
## [70] codetools_0.2-19 fastmatch_1.1-4 leiden_0.4.3.1
## [73] glue_1.7.0 downloader_0.4 ggfun_0.1.5
## [76] data.table_1.15.2 treeio_1.24.3 vctrs_0.6.5
## [79] png_0.1-8 spam_2.10-0 gtable_0.3.5
## [82] assertthat_0.2.1 cachem_1.1.0 xfun_0.43
## [85] S4Arrays_1.0.6 mime_0.12 tidygraph_1.3.1
## [88] survival_3.5-8 DElegate_1.2.1 SingleCellExperiment_1.22.0
## [91] pheatmap_1.0.12 iterators_1.0.14 fitdistrplus_1.1-11
## [94] ROCR_1.0-11 nlme_3.1-164 ggtree_3.13.0.001
## [97] bit64_4.0.5 RcppAnnoy_0.0.22 evd_2.3-6.1
## [100] GenomeInfoDb_1.36.4 bslib_0.6.2 irlba_2.3.5.1
## [103] vipor_0.4.7 KernSmooth_2.23-22 DBI_1.2.2
## [106] colorspace_2.1-0 ggrastr_1.0.2 tidyselect_1.2.1
## [109] bit_4.0.5 compiler_4.3.0 SparseM_1.81
## [112] DelayedArray_0.26.7 plotly_4.10.4 shadowtext_0.1.3
## [115] lmtest_0.9-40 digest_0.6.35 goftest_1.2-3
## [118] spatstat.utils_3.0-4 rmarkdown_2.26 XVector_0.40.0
## [121] htmltools_0.5.8 pkgconfig_2.0.3 sparseMatrixStats_1.12.2
## [124] MatrixGenerics_1.12.3 highr_0.10 fastmap_1.2.0
## [127] rlang_1.1.4 htmlwidgets_1.6.4 shiny_1.8.1
## [130] farver_2.1.2 jquerylib_0.1.4 zoo_1.8-12
## [133] jsonlite_1.8.8 BiocParallel_1.34.2 GOSemSim_2.26.1
## [136] RCurl_1.98-1.14 magrittr_2.0.3 GenomeInfoDbData_1.2.10
## [139] ggplotify_0.1.2 dotCall64_1.1-1 munsell_0.5.1
## [142] Rcpp_1.0.12 evmix_2.12 babelgene_22.9
## [145] ape_5.8 reticulate_1.35.0 truncdist_1.0-2
## [148] stringi_1.8.4 ggalluvial_0.12.5 ggraph_2.2.1
## [151] zlibbioc_1.46.0 MASS_7.3-60.0.1 plyr_1.8.9
## [154] parallel_4.3.0 listenv_0.9.1 ggrepel_0.9.5
## [157] forcats_1.0.0 deldir_2.0-4 Biostrings_2.68.1
## [160] graphlayouts_1.1.1 splines_4.3.0 tensor_1.5
## [163] locfit_1.5-9.9 igraph_2.0.3 spatstat.geom_3.2-9
## [166] cubature_2.1.0 ggsignif_0.6.4 RcppHNSW_0.6.0
## [169] evaluate_0.23 foreach_1.5.2 tweenr_2.0.3
## [172] httpuv_1.6.15 RANN_2.6.1 purrr_1.0.2
## [175] polyclip_1.10-6 future_1.33.2 scattermore_1.2
## [178] ggforce_0.4.2 broom_1.0.5 xtable_1.8-4
## [181] tidytree_0.4.6 RSpectra_0.16-1 rstatix_0.7.2
## [184] later_1.3.2 gsl_2.1-8 aplot_0.2.3
## [187] beeswarm_0.4.0 memoise_2.0.1 cluster_2.1.6
## [190] powerTCR_1.20.0 globals_0.16.3